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Abstract. We define the concept of instability index of an isolated eigenvalue 
of a non-self-adjoint operator, and prove some of its general properties. We also 
describe a stable procedure for computing this index for Schrodinger operators 
in one dimension, and apply it to the complex resonances of a typical operator 
with a dilation analytic potential. 
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1. Introduction 

In some earlier papers we showed that typical non-self-adjoint Schrodinger opera- 
tors H exhibit spectral instability in the following sense. For any e > there exist 
many A G C and / G Dom(H) such that 

||ff/-A/|| 2 <e| 



even though A is not near the spectrum of H . This behaviour occurs for the 
harmonic oscillator with a nonreal coupling constant as well as for many non- 
self-adjoint anharmonic oscillators. There is a rapidly growing literature on pseu- 
dospectral theory, which was invented to explore just such possibilities, [H, 0, 0, 0, 



ra, m m, m m m 20 . 
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In this paper we return to the same type of operator, but measure spectral insta- 
bility by a method which provides more precise information about the instability 
of individual eigenvalues. We have computed the so-called instability indices of 
the first 100 eigenvalues A n of the harmonic oscillator, and see that they appear 
to increase exponentially with n. We have carried out a similar but more limited 
exercise for the resonances of a typical Schrodinger operator with dilation analytic 
potential, and report our conclusions. 

Our definition of the instability index of an isolated eigenvalue A of if of multiplicity 
1 involves the fact that the eigenfunction / of H associated with A is different from 
the eigenfunction /* of H* associated with its eigenvalue A. In Section 2 we show 
that the instability index 

K(A) - W>T 

of A is equal to the norm of the spectral projection P associated with A, and also 
present a number of other theoretical properties of the index. 

If H := —A + V acts in L 2 (R N ) where V is a complex-valued potential then 
H* = — A + V and f* — f. Hence for any isolated eigenvalue we have 

(1) «(A) - IrN l/|2 



Ir n f z 



Note that when k(X) is large the eigenvalue is very unstable under small pertur- 
bations of the potential, and hence also unstable because of rounding errors in the 
computation. The numerical task we face is to compute the eigenvalue and the 
instability index accurately in situations in which the denominator of (1) is very 
small because the complex-valued eigenfunction / is oscillating rapidly. 

Because the spectral instability develops so rapidly as n increases, we have had to 
take great care to use computational methods which are reliable. Fortunately for 
our first problem there are independent methods of checking the values which we 
have obtained. In the second case we use the experience gained by the first problem, 
and have checked the reliability of the conclusions under the variation of several 
different parameters in the computational method. In Section 6 we summarize the 
conclusions of our investigation. 
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2. The Instability Index 

If A is an isolated point of the spectrum of a closed operator A in a Hilbert space 
H, the spectral projection P associated with A is defined by 



where 7 is any sufficiently small closed contour winding around A. The assump- 
tion that this projection has rank 1 is stronger than the assumption that A is an 
eigenvalue of multiplicity 1. 

Lemma 1. If f and f* are the normalised eigenvectors of A and A* associated 
with the eigenvalues A and A respectively, and if P has rank 1, then (f, /*) 7^ 
and the instability index of X is equal to \\P\\. 

Proof We have 

Ker(P) = (RaniP*)) 1 - 

= {<?:<<?, /*> = 0}. 

Since / ^ Ker(P) we see that (/, /*) 7^ 0. It is now easy to verify that P is given 
by 

p h -=^iif 

and hence that 

Theorem 1. If P has rank 1 then 

k(X) = sup{ ai (V) : ||V|| < 1} 

where ai(V) is defined to be the coefficient of s in the expansion of the eigenvalue 
of the perturbed operator H + sV: 

A(s) = A + ai(^)s + a 2 (V> 2 + --- 

Proof By standard arguments in perturbation theory |14| we have 

(H + sV)(f + sg + ...) = (A + sfji + ...)(/ + sg + ... ) 
where the perturbed eigenvector is normalised by 

(f+ sg +..., n = (f,n. 
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We deduce that Hg + Vf = \g + /if and (g, /*) = 0. Therefore 



(Hg,n + (vf,n=Mn 



and 



H = 



(VfJ*) 

</,/*> ' 



The proof is completed by the observation that 



sup{|<v7,r>|:||vi<i} = n/ii nr 



The instability index is also related to the notion of pseudospectrum, which is a 
geometric way of looking at resolvent norm properties. Namely if e > we put 



The sizes of these sets, which all contain the spectrum of A, measure the spectral 
instability of A. One always has 



but the RHS is often much larger than the LHS. The following theorem states that 
if k(X) is large then the component of Spec e (A) containing A is large in a related 
sense. Several similar results can be obtained in the same manner. 

Theorem 2. Suppose that the spectral projection P associated with the isolated 
eigenvalue A of A has rank 1 . Let 7 be a closed contour surrounding the connected 
component of Spec e (A) which contains A but does not intersect Spec £ (A). Then 



Spec £ {A) := {z e C : \\(z - A)' 1 ]] > a' 1 }. 



{z : dist{z, Spec(A)} < e} C Spec £ (H) 



7| > 2vre/t(A) 



where \ j\ is the length of '7. 



Proof We have 
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3. The Computational Procedures 




acting in L 2 (R). Let us first outline briefly the method for determining the eigen- 
values of H. For any z [0, oo) there exist two solutions f± of 



which vanish as x — > ±oo respectively. We introduce transfer functions g± : = 
f±/f± satisfying nonlinear first order differential equations to be given below and 
proper initial conditions at ±oo. We wish to solve the Cauchy problem for g~(x) 
for x > X_ and for g+(x) for x < X + where X_ < 0, X + > 0, and X + 

are sufficiently large. Provided we know initial conditions for g± at X± which 
correspond to f± vanishing at ±oo respectively, the two Cauchy problems can be 
solved by a standard numerical method. The question how to transfer the so-called 
admissible boundary conditions from singular points (which are ±oo in our case) 
has been extensively studied by A. A. Abramov and his collaborators (a survey of 
their results can be found in Naturally, the behaviour of the potential V has 
to be taken into account: say, if V is bounded and vanishes rapidly at infinity then 
g± are asymptotically constant as x — > ±oo, respectively. Later on in this section 
we shall consider this and other cases applying the ideas of @ . 

To locate the eigenvalues in terms of the transfer functions we choose a G R and 
consider 



This function is meromorphic on C\[0, oo) with zeros at the eigenvalues of the 
operator H. The zeros are independent of a but F may also have poles which 
depend on a. They are at the eigenvalues of the restrictions of H to L 2 (a, oo) and 
L 2 (— oo, a) subject to Dirichlet boundary conditions at a in both cases. 

To determine the zeros of F numerically we use the argument principle (cf. || for 
a contour integration procedure) to obtain their approximate positions followed by 
some variant of Newton's method to obtain accurate values. One has to be careful 
not to choose a value of a for which there is a pole close to the zero of interest, so 
it is recommended that a few different values of a are investigated. 

The numerical elaboration of the above ideas involves a substantial amount of 
preliminary work. One may find the approximate location of the eigenvalues and 
of the maxima of the eigenfunctions by an independent method. For example if one 



(3) 



Hf = zf 



F(z) = F(z; a) := g+(a) - g-{a). 
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discretises a large enough interval in the real line then one can find approximate 
eigenvalues by a standard matrix eigenvalue routine; MATLAB is ideal for this 
purpose. Another possibility is to use JWKB asymptotic formulae which also 
enable one to find an interval (X_,X + ) outside which the relevant eigenfunction 
is negligible, and a point xq at which the modulus of the eigenfunction takes its 
maximum value (see the next subsection in this regard). 

If the potential V is an even function with respect to reflection about the origin, 
then every eigenfunction is either even or odd, and the problem may be decomposed 
into two independent problems on [0,oo), with Dirichlet or Neumann boundary 
conditions at 0. This is the case in our examples. From now on in this subsection 
we concentrate on the symmetric problem and take advantage of symmetry. Actual 
computational formulae are given below for the case of the half-line. Note that 
the theory of admissible boundary conditions based upon asymptotic analysis of 
solutions of the differential equation at ±00 applies in a generic situation. 

Auxiliary Cauchy problems to be solved numerically are as follows. Let us first 
assume that V(x) ~ cx 2 , x — > 00, where c <G C is constant — this corresponds to 
the harmonic oscillator problem and its perturbations studied in the next section. 
For a fixed value of z we consider the solution /+ of (3) vanishing as x — > 00 and 
introduce a new transfer function a + satisfying 

-a' + + y/ca 2 + + ^ T (y/ca + -V+z) = 0, a+{X) 1+0 f-^) , X -> 00 

X X y C V A / 

where X is sufficiently large. According to 0, 

f' + (x)/U(x) = yfcxa + (x) 

for such x > that a + (x) exists. Moreover, one can work out the coefficients of 
the asymptotic expansion of a + for a particular potential V(x). For V(x) = cx 2 , 
say, we replace the condition at infinity by 

d y/aP-d z-y/c 

a+ix) = ~x 2 "v^' = 

choosing X so that the last term in the above formula is negligible. Along the lines 
of the mentioned paper we pose the so-called admissible condition at infinity for 
the considered potential. The above initial condition is equivalent to the boundary 
condition /+ — > at infinity, up to terms of order 0(yrg) as X — > 00. 

Next, we are going to consider the operator H with a potential vanishing rapidly 
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enough at infinity (see Section 5). Following the same approach, we introduce 

g+(x) = f + (x)/f + (x), x>0, 
and for this function obtain the singular Cauchy problem 

g' + + g 2 + -V + z = 0, g + {x)~i^, x -> oo. 

Clearly, the initial conditions vary for different types of potential. Still, for each 
particular choice of V we are able to apply the developed theory and set appropriate 
initial conditions for transfer equations. After that has been done we integrate 
those equations numerically from X from right to left to some (fixed) a > 0. The 
transfer functions a + , g + actually take their values in the Riemann sphere, so we 
have to switch between them and their inverses at certain values of x. As soon as 
\a + \ or \g + \ becomes greater than a prescribed constant we change to a^ 1 or g^ 1 , 
respectively and from this point on integrate analogous equations starting with 
proper initial conditions for the inverse functions. After a finite number of changes 
of this kind we reach the chosen point a. 

To complete the transfer procedure, consider the solution f of (3) satisfying 
/£(0) = 0. We denote 

fo( x )/fo(x) = g (x) 

and solve 

g' + g 2 -V(x)+z = 0, g (0) = 

from left to right. (Obvious changes should be made when considering an odd 
solution f Q satisfying a Dirichlet boundary condition at x — 0.) Again, we switch 
from g to g^ 1 if necessary and stop at the same point a. 

Remark that the described procedure is the simplest version of the boundary con- 
dition transfer, or pivotal condensation method we have chosen for the second 
order equation. There is an extensive literature (cf. [l|, |], 11] etc.) on the transfer 



methods where much more advanced techniques are developed. Although there 
are other possibilities, in this paper we get satisfactory results implementing the 
above version. 

Finally, we have g (corresponding to either odd or even / ) an d g+ calculated at 
the same point a. If 

F(X;a) := g+(a) - g (a) = 

then A is an eigenvalue of H. Thus, we evaluate F(z; a) for certain values of z 
as described and then find the eigenvalues of interest as zeros of F(z; a). As has 
already been mentioned, we first use the argument principle (see || for computa- 
tional formulae) to locate the eigenvalues and then apply an iterative Newton-like 
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method to obtain more precise values. After an eigenvalue has been located up 
to the required accuracy, one can compute the corresponding eigenfunction by 
recovering its values from the transfer functions g + , go- 

The method proposed has been implemented as a universal Fortran 77 code in- 
cluding all the basic procedures described above in this section. Auxiliary Cauchy 
problems have been solved by a standard routine based on the Runge-Kutta- 
Merson fourth order method. In our computations we have used 32-bit and 64-bit 
arithmetic. 

The question remains how to choose a — although the zeros of F(z; a) do not 
depend on a, in practical computations the choice of a does play a significant role. 
In fact, we investigated different functions F(z; a) for a wide range of a. One of the 
possible choices is a = 0. For problems with even potentials the zeros and the poles 
of F(z; a) interchange (see Figures la, lb, 2a). We use contour map plots to find 
initial guesses for eigenvalues when there is no other a-priori information about 
their location. Plotting contour maps with the use of Matlab 5.2 has also helped us 
to avoid poles when looking for zeros. Thus, in a generic situation we recommend 
that one calculates the values of F(z) (for which we have used our Fortran code), 
then produces the plots (Matlab graphics) and, finally, finds eigenvalues accurately 
(a standard iterative Fortran procedure). 

We regard our numerical results obtained via the above method as reliable. In 
particular, this is confirmed by several values of a providing entirely different func- 
tions F(z) whose zeros coincide. These results are to be reported below in the 
following two sections. 

3.2. Calculating the Instability Index. The second stage in the process is 
to compute the instability index defined by (1). The obvious method, namely 
calculating the two integrals after first determining the eigenfunction numerically, 
is highly inaccurate if the instability index is large. The reason is that the integrand 
in the denominator is highly oscillatory, and the evaluation of such integrals is 
problematical. The following method is much superior in applications. Below we 
present a technique suitable for an arbitrary (not necessarily symmetric) potential. 

We introduce four functions as follows. For x G [a, X + ] we define 
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where / is the eigenfunction associated with the eigenvalue A. Similarly for x G 
[X_ , a] we define 



It is obvious that 



h-(x) := f(x)- 2 jT f(s) 2 ds 
k.(x) := \f(x)\- 2 T \m\ 2 ds. 

J .A — 

fc_(a) + fc+(a) /l + |/(x)| 2 rfx 



|/i_(a) + /i+(a)| f*+f(x) 2 dx 

which converges exponentially rapidly to /t as X + — > +00 and X_ — > — 00. The 
task is to find a procedure to evaluate the four functions accurately. We consider 
only h-i the others being similar. It follows from its definition that h-(X-) = 
and that h- satisfies the differential equation 

(4) h-(x)' = 1 -2g_(x)h_(x). 

This may be solved numerically, say, by a Runge-Kutta method to determine h- (a). 

It is important to be sure that the solutions of (4) and the other three equations are 
stable. It suffices to note that Re g + (X + ) < and Re g_(X_) > 0, which implies 
the stability of the solutions h + , k + and h-, k_ from right to left and from left to 
right, respectively. This has been confirmed by numerical tests. The same is true 
for the transfer equations quoted in the previous subsection — the solution a + , 
for instance, is known to be stable from right to left which is essential for practical 
computations. 

There is a potential problem in that if f(b) = for some b G (X_,a) then h_ is 
usually infinite at that point. Generically one does not expect a complex- valued C 2 
function of a real variable to vanish anywhere, and we have not seen this problem 
arise, but one needs to discuss how the method should be adapted in the event 
of its occurrence. There are two cases, which are distinguished numerically by 
whether h-(x) — > as x — > b or — > 00 as x — > b. Note that since / is a 
non-zero solution of (3) and f(b) = it follows that f'(b) 7^ 0, so |<7(x)| — > 00 as 
x — > b. 

Lemma 2. If f(b) = and f h x _f{s) 2 ds = then h_(x) -> as x — > b and 
g(x)h-(x) — > 1/3 as x — > 6. 
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Proof Neglecting lower order terms we have 

h-(x) ~ (x — b)/3 
g(x) ~ (x-by 1 

as x — > b. The results follow. 

Thus, in this case we are still able to integrate the same equation (4); the point b 
is, in fact, regular rather than singular. 

The more standard case is that in which f(b) = and /£_ f(s) 2 ds ^ 0. Clearly 
— > oo as x — > 6. 

Lemma 3. // w;e put h-(x) := £/ien h-(x) — > and h-(x)g(x) — > as 

x — > 6 ana 1 

(5) /L(ar)' = /i-(x) 2 -2/i_(x)^(x) 

/or a// x near 6. 

Proof Neglecting lower order terms we have 

/W(z-&) 2 



fe_(x)^(a;) 



f_W{x-b) 
fx. f(s) 2 ds 



as x — > 6. The verification that h-(x) satisfies the differential equation (5) is 
routine. 

Thus, in the considered case we recommend to change to h at x — b and integrate 
(5) instead of (4). 

Naturally, the stability of the procedure proposed in this subsection depends heav- 
ily on a (though the exact value of k does not depend on the norm of /). A proper 
choice of a is very important and can essentially influence the results. Choosing 
a = argmax\f(x)\ seems to be a reasonable way. 

Compared to standard approaches the above mentioned technique has two clear 
advantages. First, we do not need to evaluate the fast oscillating integrands f(x) 2 
and |/(a;)| 2 themselves — instead, we integrate several auxiliary ODEs. Secondly, 
this procedure is numerically stable. In the cases which we have examined the 
solutions h±, k± change quite slowly and smoothly. 
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3.3. Possible Difficulties. If the instability index of an eigenvalue is very large 
then it is clear from Theorem 2 that the eigenvalue is intrinsically difficult to 
compute. One mechanism by which this can occur in computations is that at the 
eigenvalues, for which one knows that F(X) = 0, one also finds that F'(\) is very 
small, so it is not possible to locate A accurately. The following theorem provides 
a link in one direction between these two phenomena at a theoretical level. 

We assume that 



on [a, oo) subject to g+(z, x) ~ iy/z as x — > oo. Let g~(z, x) be the solution of the 
same equation on (— oo,a] subject to g-(z,x) ~ —iy/z as x — > — oo. We put 



as usual so that F(X) = if and only if A is an eigenvalue of H. 

Theorem 3. Let k(X) be the instability index at an eigenvalue X, let f be the 
associated eigenfunction and assume that g := f'/f is bounded on R. Then 




oo. Given a G R and 



z G C satisfying Re (iy/z) < 0, let g + (z,x) be the solution of 

g'(x) + g(x) 2 + z- V(x) = 



F ( z ) = 9+(z,a) -g-(z,a) 



K (A)|F'(A)|||s||oc>l. 



Proof If e > is small enough there exists 



// = A + 



e 



+ 0(e 2 ) 



F'(A) 



such that F(fj) = e. Now put 



g + (x) := 
for the appropriate values of x, so that 



g-{x 



g-(n,x)+e/2 
g+(fi,x) -e/2 



9+(a) ~9-( a ) = °- 



We have 




V(x) - ii- g-(x) 2 
V(x) -fi-g + {x) 2 
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under the following conditions on V. If x > a then 

V(x)-V{x) = g' + (x) +/2 + g + (x) 2 - V(x) 

= g + (ii,x)' + 11+ (g+(n,x) - e/2) 2 -V(x) 
= ~eg+(fi,x)+e 2 /4: 

while if x < a we must have 

V(x)-V(x) = gi(x)+i2 + g_(x) 2 -V(x) 

= g-(n,x)' + fi + (g.(fi,x)+e/2) 2 -V(x) 
= eg-(fj.,x) +e 2 /4 

Therefore 

\\V - VWoo = sup{e\\g + (ii, -)||oo + 0(e 2 ),e\\g.(fi, -)||oo + 0(e 2 )}. 

But g±(/i,x) — > g(x) uniformly as /i — > A by the assumptions of this section, so 

||V"- V||oo =e||^||oo + o(e). 

The statement of the theorem now follows from the formula for the instability 
index given in Theorem 2. 

In the two examples considered below, F'(z) is very small for large values of \z\, 
so it is impossible to determine its zeros. This seems to be the main barrier to the 
determination of large eigenvalues. 

4. The Harmonic Oscillator 

4.1. Basic Facts. Consider the operator H defined by (2) with the potential 
V(x) = cx 2 , referred to as H Q in the rest of the paper. The eigenvalue problem 
for H is called the harmonic oscillator problem and is known to have infinitely 
many eigenvalues \^ = y/c(2n + 1), n = 0,1,.... The corresponding eigen- 
functions f n = C n e~^ x l 2 4> n {\/~cx) where C n are normalising constants, 4> n denote 
Hermite polynomials, Re y/c > 0. These eigenfunctions are either even or odd: 
fzk(x) = f2k{~x), f 2 k+i(x) = -f 2 k+i(-x), k = 0, 1, ... . As proposed in Section 3, 
we consider H on the half-line adding either Neumann or Dirichlet boundary con- 
ditions at the origin. We have used it as a sample problem to check the above 
method for finding eigenvalues and eigenfunctions. Indeed, the results thus ob- 
tained are in very good accordance with the theory; they confirm the reliability of 
the method. It allowed us to calculate ~ 100 eigenvalues of H Q up to the accuracy 
5 G (10~ 10 , 10~ 4 ). 
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When implementing the method of Section 3 we found that the accurate numer- 
ical determination of the eigenvalue A n for n > 100 is not possible using double 
precision (64-bit) arithmetic, particularly because F'(z) can be very small near 
the points where F(z) = 0. We computed the instability indices for the first 100 
eigenvalues, using the JWKB approximation to the eigenfunction f n associated 
with the eigenvalue A^ as described in || [TIJ . This approximation suggests that 
takes its maximum near x = a n and 

where the real constants a n and r\ n are computed from 

v^(2n + l)=^ + m 2. 

Having an appropriate value of a n is, of course, helpful when we calculate K n by 
means of the method given in Subsection 3.2. 

4.2. Perturbations of the Operator H a . Let us present some results concerned 
a perturbation of the harmonic oscillator operator 

H w :=H + W(x). 

We have investigated various perturbations of the form W(x) = ee tmx for a range 
of fairly small e. The reason is that looking for the perturbation providing the 
most unstable results, one has to choose W(x) as follows. From the perturbation 
theory formula cited in the proof of Theorem 1 one can easily see that among 
perturbations satisfying 

|W(a;)|<e 

the function 

provides the worst perturbation of the n-th eigenvalue of H Q . Indeed we then have 
(6) X n = X^+eK(X^) + o(e). 

If we only take account of the first term of the JWKB expansion for f n , we obtain 
the perturbing potential 

W n (x) = ee 2iVnX 



after removing an irrelevant phase factor. The expectation that W n provides a 
perturbation of the eigenvalue almost as great as that due to W n is tested below. 
We tabulate below the absolute values of the corrections to several eigenvalues 
A n of H w calculated numerically by means of the method described in Section 3. 
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Tables 1-3 contain the values of - A n | for c = \fi and W(x) = ee imx . The 
figures related to e = give the absolute errors of the computation of the n-th 
eigenvalue of H a . 



Table 1. Values of 
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e\m 


1.0 


5.0 


6.0 


6.4133 


7.0 


10.0 





io- io 












io- e 


7.9 • IO" 7 


6.6 • 10~ e 


7.7- IO" 6 


8.0 • 10~ e 


7.8 • 10" e 


2.9 • IO -7 


10~ 5 


8.6- 10-" 


6.4 • IO -5 


7.9 • IO -5 


8.1 • IO" 5 


7.7 • IO -5 


1.9 • 10~ e 


10~ 4 


8.7- IO" 5 


6.4- IO" 4 


8.0- IO" 4 


8.0 • IO" 4 


7.7 • IO -5 


1.8 • 10~ 5 


io- 3 


8.7- 10~ 4 


6.42 • 10~ 3 


7.97 • 10~ 3 


8.16 • 10~ 3 


7.66 • 10~ 3 


1.7- IO" 4 



Table 2. Values of \X$ - X n \,n = 19, 2r] n = 9.1884 



e\m 


5.0 


9.0 


9.1884 


10.0 


20.0 





io- 8 










IO" 8 


1.0- IO" 7 


4.0 • l(T e 


3.9- 10- s 


1.1 • io~" 


6.1 • IO" 7 


IO" 7 


2.6- 10- e 


3.25- 10~ 5 


3.22- 10~ 5 


2.63- 10~ 5 


6.2 • IO" 7 


io- 6 


3.7- 10~ 5 


3.25- IO" 4 


3.20- IO" 4 


2.77- 10~ 4 


6.1 • IO" 7 


IO -5 


4.0- IO" 4 


3.20 • 10~ 3 


3.20 • IO -3 


2.80 • 10~ 3 


6.3 • IO - '' 



Table 3. Values of \X^ - X n \,n = 29, 2i] n = 11.3014 



e\m 


10.0 


11.0 


11.3014 


12.0 


20.0 





io-" 










io- io 


2 -10-" 


2 • 10- e 


3- 10- e 


3 ■ 10- e 


io- B 


IO" 9 


1.2 • IO -5 


1.5 • 10" b 


1.4 • IO -5 


1.4 • 10" b 


io~ e 


IO" 8 


1.14- 10~ 4 


1.48 • 10~ 4 


1.46 • IO" 4 


1.34- IO -4 


io~ B 


io- v 


1.141 • 10~ 3 


1.454- IO -3 


1.441 ■ 10~ 3 


1.336 ■ IO -3 


io~ e 


10~ e 


0.01142 


0.01455 


0.01453 


0.01334 


10~ d 


IO" 5 


0.11657 


0.14927 


0.14563 


0.13670 


10~ d 



First of all, the above results show that the values of \X$ — A n | are approximately 
proportional to e, that is, confirm formula (6) numerically. In fact, for very small 
e the method can only feel the first order corrections to the eigenvalues of the 
harmonic oscillator within the chosen accuracy as expected. Secondly, maximal 
perturbations of eigenvalues are observed for m ~ 2i] n which justifies the above 
arguments. 
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4.3. Another Approach. Instability Index. We have also investigated the 
harmonic oscillator using the quantum mechanical creation and annihilation oper- 
ators A* and A. This is not possible for generic differential operators, but provides 
a method of testing the general algorithms developed in the last section. In this 
language 

H = P 2 + cQ 2 

= -(A* - A) 2 /2 + c(A* + A) 2 /2 

= (c - l)A* 2 /2 + (c + I) A* A + (c - l)A 2 /2 + (c + l)/2. 

If {0 n }^Lo i s the orthonormal basis of Hermite functions in H := L 2 (R), then 
A(p n = y/n-4>n-i and A*(f) n = \Jn + 10 n +i f° r an n i an d we may represent H by 
means of the infinite matrix 



H ,m,n ■ < 



a m if m = n 

b m if n = m + 2 

6 n if m = n + 2 

otherwise 



with respect to this basis, where m,n — 0,1,2, .. . and 

a m := (c+l)(m + l/2) 

b m := (c-l){(m+l)(m + 2)} 1 / 2 /2. 

The even and odd subspaces 7io and 7Yi with respect to reflection about are 
invariant under H Q , and these subspaces may also be characterised by 

H = lin{4>2n ■ n = 0, 1, . . . } 
Hi = lin{<p2n+i ■ n = 0, 1, . . . }. 

Restricting the matrix to either of these subspaces renders it tri-diagonal, so numer- 
ical computations are particularly easy and accurate. We compute the instability 
index of an eigenvalue A r = ^/c(2r +1) for r = 0, 1, . . . by evaluating 



Kir . — 



where / is the eigenvector associated with A r , obtained by solving the obvious 
recurrence relation starting from n — 0. Note that A r is taken to be an exact 
eigenvalue of the infinite matrix, not an eigenvalue of the truncated N x N matrix. 
For a particular eigenvalue X r , N must be large enough for the coefficients /„ 
with n > N to be insignificant, but not so large that the recurrence relation 
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becomes unstable. For r > 50 it is not possible to satisfy both of these conditions 
simultaneously using standard double precision 32-bit arithmetic, and we used the 
high precision arithmetic of Maple V.4. 

The delicacy of the computations is indicated by the evaluation of Kioo ~ 2.5594 x 
10 16 for c = \fi. For N = 200 this required us to use the command ""Digits = 
30' in Maple V.4, but putting N = 500, we only obtained the same result for 
'Digits = 110' or greater. The instability in the solution of the recurrence relation 
is evidently more important than the contributions of the terms of the series in the 
range 200 < n < 500. The following results (see Table 4) were all obtained with 
N = 200 and 'Digits = 100', and appear to be reliable. 

The instability indices tabulated below have been obtained in two independent 
ways. The methods developed in this and the previous sections turned out to 
provide very close results for the first 40 eigenvalues. This can be seen from Table 4 
where is related to the method of this section, and to that of Section 3. 
The figures obtained for n > 40 are clearly different for the two methods although 
they are qualitatively of the same order. 



Table 4. Instability indices of H Q) c — \fi 



n 





10 


20 


30 


40 


50 


™n 


1.0404 


14.2777 


563.2146 


2.5789-10 4 


1.2625-10 6 


6.3627-10 7 




1.0404 


14.2777 


563.2146 


2.5789-10 4 


1.2625-10 6 


6.3649-10 7 


n 


60 


70 


80 


90 


100 




™n 


3.2734-10 9 


1.708M0 11 


9.0059-10 12 


4.7860-10 14 


2.5594-10 16 






3.2922-10 9 


1.7110-10 11 


8.9063-10 12 


4.0052-10 14 


1.926M0 16 





The growth of the instability index corresponds to the values of | A^ — A n | increasing 
with n (see also Tables 1-3). Results to be cited below provide another numerical 
evidence of this fact. In Table 5 the values of — A n |, c = y/i, corresponding to 
the perturbing potentials W n are given. Comparing Table 5 to Table 4 we conclude 
that fi n ~ K n which indicates reasonably good agreement of our numerical results 
and perturbation theory. 



Table 5. Values of |A^ - A n |, W n (x) = ee 2ir >™ x , c=V~i 
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e\n 


on 

30 


a a 

40 


r - a 

50 


60 


1 A — ("> 

10 


r\ n et 1 r a o 

0.021542 


1.19860 






10 


0.002105 


0.10/4/ 






io- y 


0.000216 


0.01056 


0.54950 




lO" 9 


2.2 • 10~ 5 


0.00105 


0.05518 


2.4921 


10 -io 


3 • 10- e 


0.00010 


0.00587 


0.2587 


io- n 




10" b 


0.00059 


0.0279 


io- 12 






6 ■ 10" 5 


0.0028 





io- s 


5 • l(T e 


10" 5 


10" 4 



Analysing the rate of divergence of the instability index of H Q in Table 4, one can 
notice that it grows exponentially: K n ~ e 0An for the studied range of n. 

The eigenfunctions for the harmonic operator with nonreal coupling constant do 
not form an unconditional basis, 0. If they formed a conditional basis the projec- 
tions P n associated with the eigenvalues A n = *Jc(2n + 1) as in Lemma 1 would be 
uniformly bounded in norm by a standard argument, ||12|| . However, we have ob- 



tained strong numerical evidence that the norms increase exponentially with n. We 
therefore make the conjecture that for nonreal coupling constant the eigenfunctions 
of the harmonic oscillator do not form a conditional basis. 

More precisely let iV > and let Pn be the spectral projection of H Q associated 
with the first iV complex eigenvalues where they are ordered in increasing 
absolute values. Explicitly 

JV-l if f*\ 

r> f \ " \J i J nl f 

^Nj ■- 2^ 77 777 Jn ' 

n=0 V"> Jn/ 

If \\PnI — f\\ — ► as N — > +oo for all / £ TC then the uniform boundedness 
theorem implies that there exists a constant C such that \\Pn\\ < C for all N. 
From the inequality 

\\Pn — Pn-i II < 2(7 

we are then able to deduce that the instability index bounded function of n. 

This conflicts with the numerical evidence that these indices increase exponentially 
with n. We have attempted to confirm the exponential increase by using the JWKB 
approximations to the eigenfunctions constructed in 0, |10[j , but the eigenfunctions 
oscillate so rapidly for high eigenvalues that the JWKB approximations were not 
useful. While we have not proved the exponential increase of K n the corresponding 
result for the pseudospectrum (resolvent norms) has been proved in [|1(J not just 
for the harmonic oscillator but for a wide range of anharmonic oscillators. 
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5. Complex Resonances 



Hf(x) :=-- r L + V(x)f(x) 



5.1. Definitions. Let H be the Schrodinger operator 

H 

dx 1 

acting in L 2 ([0, oo)) subject to Dirichlet or Neumann boundary conditions at x = 0, 
where the potential V is bounded and vanishes at infinity. For any positive constant 
c we define 



HJ{x) := (D c HD c -xf 



(7) = - C -*£L + V(cx)f(x) 

where D c is the unitary dilation operator 

D c f(x) := ^f(cx). 

We observe that H c is unitarily equivalent to H. If V is an entire function on C 
then the formula ([I]) defines a family of non-self-adjoint operators parametrised 
by c G C, c 7^ 0. Under suitable conditions the eigenvalues of these operators 
are known to be independent of c, and are called resonances of H; see [0, [L3| for 



expositions of the theory of dilation analytic resonances. Since the operators H c 
are unitarily equivalent for values of c with the same argument, we only consider 
c of the form c := e %e l 2 where < 9 < n/2. 

We investigate the particular case of the operator 



H f(x) : =-0+^ 2/fe V(x 



where b > is to be fixed. If one imposes a Dirichlet boundary condition at x = 0, 
this operator determines the evolution in the zero angular momentum sector of 
a three-dimensional quantum particle trapped by a rotationally invariant barrier, 
where the particle may tunnel through the barrier and escape to infinity. Because 
the potential is non-negative and vanishes rapidly at infinity, Hq has absolutely 
continuous spectrum [0, oo) and no eigenvalues. A direct calculation shows that 



IfnfU) : H,,, ,./V) -e-'^+e'Ve-^^/^i 

ax z 



We consider Hg subject to either Dirichlet or Neumann boundary conditions at 
x = 0. The potential of this operator vanishes rapidly as x —>■ oo provided < 9 < 
7r/2. Under this condition Hq has essential spectrum e~ 10 x [0, oo) and also some 
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isolated eigenvalues in the sector {z : — 9 < aigz < 0}, these being independent of 
9. 

For large values of b (we take b = 100) the potential of H e is similar to that of 
the complex harmonic oscillator, and the eigenvalues of Hq are close to the values 
{2n +1: n — 0, 1, . . . }. For smaller values of b there are several resonances very 
close to the positive real axis, but at a certain point they turn sharply away into 
the lower half plane. 



5.2. Location of Resonances. The reason for there being resonances very close 
to the real axis is as follows. Let us consider the operator Hq as a perturbation of 
the harmonic oscillator operator H Q . In our notation we now have 

H = H O + W(t, v), W(t, u) = t 2 (e~ t2 / h ' 2 - 1) 

where we put t = e ld / 2 x and v = 1/b 2 . Regarding v as a small parameter we 
expand 

w(t, u) = y: w k (ty = £ { —^ } — v k . 

k=l k=l K - 

Again, for an arbitrary n, following the standard perturbation theory approach, 
we expand the n-th eigenvalue of H e as 



(8) A»~Ai°) + ]>>^ 



A" 

k=l 



which is a non-convergent asymptotic expansion, and calculate 
(9) „ = l^Jp = -CI f ,e-Ulm- 

J-oo Jn al J -°° 

Here we follow the notations of Section 4: /„ are the eigenfunctions of H Q and 4> n 
are Hermite polynomials. 

Formula (9) implies that the first order correction is real and does not depend on 
9. The same is true for all fik- Indeed, it is easily seen that the parameter 9 enters 
the problem in a specific way. If one passes to the new variable t and proceeds 
with calculation of higher order corrections, all the relations thus obtained do not 
contain any complex values except for t as an integration variable. Thus, one only 
deals with integrals of the form Jf^ p(t)dt which do not depend on 9 and, therefore, 
are real. 
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Using the creation-annihilation technique based on the corresponding decompo- 
sition of the operator H (see the previous section), we calculate the first order 
correction for the n-th eigenvalue implicitly. Thus, formula (9) becomes 

(10) in =m(n) = -|j(2n 2 + 2n + l), n = 0, 1, . . . . 

Following the same numerical procedure (see Section 3) we compute some of the 
eigenvalues of H e . These results can be found in Subsection 5.4. 

5.3. Numerical Range and Complex Resonances. The resonances must turn 
away from the real axis as their absolute value increases, because of the fact that 
a resonance z is an eigenvalue of He- This implies that 

0<6»<tt/2 

where N(9) is the numerical range of the operator Hq. The numerical range of He 
is defined by 

N e := {(H e f,f): ||/|| = 1} 

= {J{V (x)\f{x)\ 2 + e^\f(x)\>}dx: ||/|| = 1} 

C {jv e (x)\f(x)\ 2 dx: \\f\\ = 1} + e -*[0,oo) 
C conv{Vg(x) : x G R} + e~ ie [0, oo), 

where 

V s (x) :=e id x 2 e- eZ6x2 l h2 . 

For small positive 9 the set N(9) crosses the real axis near x ~ 0.4616 2 . 

Indeed, if we represent V e (x) = X e (x) + Y e (x) and denote the point where N(9) 
meets the real axis by (Bg, 0) then a simple calculation gives us 

B e = X e + Y e cot 9 = b 2 max{x 2 e" a;2 (2 - x 2 ) : x e R} + 0(9) = 0.4616 2 + 0(9). 

Therefore the imaginary parts of any resonances must start decreasing before their 
real parts reach this value. This is in good accordance with the numerical data 
quoted in the next subsection. 
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5.4. Numerical Results. Lower eigenvalues of the operator H g lying close to 
the real axis for different values of v are given in Tables 6, 7. The computed 
eigenvalues proved not to depend on 9, so our numerical results are in agreement 
with the theoretical arguments of Section 5.1. The fact that for a range of 9 lower 
eigenvalues coincide up to a high accuracy shows the stability of our method as a 
whole. 

Remark that the results to be reported below are consistent with formulae (9) and 
(10); they confirm, in particular, that ji\ < 0. On the other hand, these results 
illustrate the fact that series (8) is asymptotic rather than convergent. This only 
implies that the imaginary parts of the resonances have to be very small within 
the regime for which the asymptotic expansion provides useful information. 



Table 6. Resonances of H , n = 0, 1 



V 


Ac 


Ai 


0. 


1. 


3. 


10" 4 


0.999925 


2.999677 


4 • 10" 4 


0.999700 


2.998502 


10" 3 


0.999251 


2.996253 


io- a 


0.992475 


2.962115 


0.04 


0.969405 


2.824312 


0.1 


0.920295 - 2. • 10 _B i 


2.560861 - 0.003347i 


0.2 


0.822647 - 0.005282i 


2.028250 - 0.249944i 


0.25 


0.768023 - 0.019417i 


1.850388 - 0.425748i 



Table 7. Resonances of H , v = 10 4 
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n 




n 


\ 

An 


O 




0.999925 


1 O 

18 


' ) / • A / ior'7i 

36.948571 


1 


;\ ooo/^TT 

2.999077 


OA 

20 


40.936844 


o 

2 


4.999025 


no 

22 


/i /i oo ooo r 

44.923925 


3 


f* aaoi or 

6.998125 


24 


A O AAA^AT 

48.909797 


A 

4 


o oor*oo /i 

8.996924 


or* 

26 


ro oo/i/ir*o 

52.894462 


5 


1 O OO CT O 

10.995877 


28 


C/' 1 CWOOO 

56.877922 


6 


1 o ru v ) < • o 1 ) 

12.993623 


QO 

30 


60.860175 


( 


1 /i on 1 o /i /i 
14. yy 1344 


4U 


oU. /533ZI 


8 


16.989119 


50 


100.616298 


9 


18.986242 


60 


120.449097 


10 


20.983418 


70 


140.237430 


12 


24.976502 


80 


160.000434 


14 


28.968399 


90 


179.874965 


16 


32.959086 


100 


199.664121 



Along with Table 7 we present some plots (see Figures la-lc). They include 
contour maps of the function F(z; a) defined in Section 3 whose zeros are the 
eigenvalues we are looking for. One can see that the zeros and the poles of F 
interchange (we have plotted F for Neumann boundary condition at x — 0, i.e., its 
zeros are the even eigenvalues, while the poles correspond to the odd ones). For 
\z\ < 200 we have discovered 100 eigenvalues all being real up to the accuracy 
5 = 10 -6 . Note that for different values of a (an intermediate matching point) 
we obviously get quite different functions F(z; a) (compare Figure la to lc) while 
their zeros remain the same. 

Given a certain number n we watch \ n changing as v increases and compare this 
eigenvalue with its first order approximation (see (8)). 



Table 8. Eigenvalues of H d , n = 10, 30 



24 A. ASLANYAN AND E. B. DAVIES 



V 




\{p) | /m\,, 




\{P) /QA\, 

A 30 +^i(30)z/ 


1 n— 4 

10 


20.9834 


on noo/i 

20.9834 


60.8602 


60.8604 


1 n-3 

1U 


on oQon 
20.8332 


on oo a o 
20.8343 


o9.o / 81 


59.0043 


2 • iir 3 


20.6643 


20.6685 


58.0976 


58.2085 


4 • io~ 3 


20.3195 


20.3370 


54.9090 


55.4170 


6 • KT 3 


19.9644 


20.0005 


51.2492 


52.6255 


8 • KT 3 


19.5531 


19.6740 


46.2597- 0.1492i 


49.8340 


9 • 1(T 3 


19.3747 


19.5083 


43.3597- 1.5496i 


48.4383 


io- 2 


19.2045 


19.3425 


41.2601 - 3.2488i 


47.0425 



It is seen from Table 8 that only for a narrow range of v do the perturbation theory 
formulae (8) approximate actual eigenvalues (compare to Table 6). As v increases 
a typical eigenvalue deviates gradually from the value given by (8) and at some 
stage its imaginary part becomes substantial. 

The values of the instability indices of resonances depend on 9, even though the 
positions of the resonances do not. We have observed that the indices are in fact 
monotonically increasing functions of 9. While this is not surprising we have no 
proof of the fact. We have also observed that the instability indices are increasing 
functions of b, provided one follows the 'same' resonance as b increases. The insta- 
bility indices K n , n = 0, 10, 20, computed for a wide range of 9 and b are given in 
the following three tables. 
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Table 9.1. Instability indices K n , n = 



25 



v \e 


tt/30 


tt/8 


tt/6 


tt/4 


tt/3 


tt/2.5 


0. 


1.002750 


1.040381 


1.074570 


1.189207 


1.414214 


1.798908 


10~ 4 


1.002750 


1.040378 


1.074563 


1.189185 


1.414134 


1.798588 


io~ 2 


1.002729 


1.040039 


1.073886 


1.186951 


1.406329 


1.769120 



Table 9.2. Instability indices K n , n = 10 



v \e 


tt/3000 


tt/300 


tt/30 


tt/20 


tt/10 


0. 


1.000028 


1.0233 


1.3299 


1.8249 


6.6784 


10~ 4 


1.000027 


1.0209 


1.3294 


1.8243 


6.6728 


io- 2 


1.000031 


1.0028 


1.3046 


1.7562 


5.9928 


v \e 


tt/8 


tt/6 


tt/5 


tt/4 


tt/3 


0. 


14.2777 


57.4539 


195.9499 


1565.2614 


1.3645551 • 10 b 


io- 4 


14.2836 


57.3505 


195.4619 


1558.9429 


1.3507237- 10 b 


io- 2 


12.3306 


45.7594 


143.8155 


965.0957 


4.6122845 • IO 4 



Table 9.3. Instability indices K n , n = 20 



v \e 


tt/20 


tt/10 


tt/6 


tt/5 


tt/4 


tt/3 


0. 


5.9275 


113.5766 


9850.7214 


1.1753- 10 b 


7.4538 • 10 s 


3.6609- 10 1U 


io- 4 


5.9190 


113.1898 


9782.2812 


1.1641 • IO 5 


7.3412 • 10 s 


3.6119- 10 1U 


IO" 3 


5.6055 


109.2017 


9128.5324 


1.0610- IO 5 


6.3605 ■ IO 6 


3.4647- 10 iU 


io- 2 


4.7349 


75.0929 


4811.8348 


4.5459 • IO 4 


1.6894 • IO 6 


1.7152 • 10 y 



Finally, let us cite some results obtained for b = 10. In this case we have found 
numerically several resonances, which are real up to the chosen accuracy 5 = 
10~ 4 , and a series of complex ones. As is seen, starting from about n = 20 their 
imaginary part rapidly increases in absolute value. In fact, for different values of 
8 the number of resonances with negative imaginary parts varies. The resonances 
and the relevant instability indices are tabulated below. In Table 10 we cite the 
eigenvalues A n of Hg along with the corresponding instability indices k„ calculated 
for 6 = tt/4 and 6 = tt/16. 

The data of Table 10 is illustrated by the plot of F(z) (Figure 2a) and its contour 
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maps (Figures 2b, 2c). Remark that the largest instability indices for 9 = n/A 
and 9 = n/16 correspond to the 26-th and the 24-th eigenvalue respectively (here 
we have concentrated on even eigenvalues; odd eigenvalues behave similarly). 
Figures 2b, 2c also show that the most unstable eigenvalues are related to n = 24 
and n = 26. They appear to be the first eigenvalues with negative imaginary parts 
- as one can see, the following eigenvalues go to the complex plane quite abruptly. 
We do not have any theoretical explanation of this fact except for the remark on 
the boundedness of n n made in the end of Section 4. Anyway, the contour maps 
and the values of n n agree very well and imply the same — the maximum of n n is 
obtained for the 'critical' range of the spectral parameter where eigenvalues start 
moving away from the real axis. 

Note that though for 9 = n/16 and 9 = ir/4 the contour map plots are quite similar, 
this only means that the first 28 eigenvalues coincide for the two operators. As 
we know, there are no eigenvalues of Hg below the line e~ %e x [0, oo). The spots 
indicating the zeros of the function F(z) which are beyond the range {z : —9 < 
argz < 0} correspond to solutions of Hgf = zf growing at infinity rather than 
decaying. Thus, in the considered example we should only regard the first 28 zeros 
as the eigenvalues of Hg, 9 = tt/16. They coincide with those obtained for 9 = tt/4 
as we expected. We believe that our results are reliable because of their stability 
under the variation of several parameters involved in the problem. 

Table 10. Values of A„ and n„ for v = 10~ 2 
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n 


\ 


n i a 
K n , = TV/4: 


f\ 1 1 P 

K n , V = 7T/16 





0.9925 


1 1 O TA 

1.1870 


1 A A AT 

1.0097 


2 


A a a a a 

4.9009 


O O A O O 

2.8983 


1 (~\P O A 

1.0684 


4 


O POOP 

8.6836 


11.3609 


1 1 AA 

l.zlUU 


6 


1 o o o r a 

12.3350 


A A A TTO 

49.4772 


1 A A ( ' O 

1.4462 


8 


15.8488 


O 1 A A 1 O A 

219.4180 


1 TA T /I 

1.7974 


1 A 

10 


19.2174 


A/" A C~ A CTO 

960.5058 


O OA 1 O 

2.2918 


1 o 

12 


OO /I O 1 o 

22.4312 


/in^r oo 

4075.82 


o c\p r o 

2.9652 


1 A 

14 


O C /I 70 


1 p c 1 c i a4 


o or fyp 

3.8576 


16 


OO O A OO 

28.3422 


/"* OO/^A 1 a4 

6.2860 • 10 


r A A o o 

5.0033 


18 


o 1 Ann /I 

31.0004 


O 1 AOA 1 A& 


A AOA 

6.4039 


20 


OO 1 O A AAAO ' 

33.7512 — 0.0003^ 


1 OATO 1 At) 

1.3978 ■ 10 


9.5706 


22 


OtX tXAAO A AA1 /(»' 

35.5098 — 0.0014? 


1 crr'rn 1 At) 

1.5650 ■ 10 


1 A A A 1 O 

10.0018 


24 


3 / .0693 — 0.1593? 


n yj roc 1 nt> 

2.4535 • 10 


10. 2o3 ( 


26 


38.7468 - 1.0004? 


2.8963- 10 e 


8.8755 


28 


39.8045 - 2.0367i 


2.5627- 10 e 


7.0743 


30 


41.2601 - 3.2488i 


1.8339- 10 e 




32 


42.6021 -4.7565i 


1.3337- 10 e 




34 


45.6102 - 8.0230i 


4.2730 • 10 5 




36 


47.0034 - 9.8515i 


2.3134- 10 b 





6. Conclusions 

The instability index of an eigenvalue of a non-self-adjoint ordinary differential 
operator was defined in Section 2, where we investigated its theoretical proper- 
ties. We have described a known general numerical procedure for computing the 
eigenvalues of the differential operator and have introduced a new and numerically 
stable procedure for computing the instability indices. 

In order to test this procedure, we have carried out extensive computations for 
the harmonic oscillator with a complex coefficient. The eigenvalues \ n of this 
operator are given by an exact formula, and we found close agreement between 
the formula and our numerical results for n < 100. We have also computed the 
instability index by two independent methods, the first being the general procedure 
mentioned above. The second uses a special numerical technique only available for 
the harmonic operator, but capable of yielding extreme accuracy if implemented 
in Maple with high precision arithmetic. The instability indices of the first 40 
eigenvalues obtained by the two methods were found to be in close agreement; see 
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Table 4. 

The discrepancies between the two methods are partly explained by the very high 
values of the instability indices of the eigenvalues A n for n > 40. This phenomenon 
was first observed for the harmonic oscillator in |J [K| where we approached the 
phenomenon via pseudospectral theory. Our current approach has the advantage 
that it provides a quantitative measure of the instability of individual eigenval- 
ues under small perturbations of the potential. We have carried out numerical 
experiments and confirmed that the size of the effects predicted matches what we 
have observed for a particular perturbation. In Section 4 we have conjectured on 
the basis of the numerical results that the eigenfunctions of the complex harmonic 
oscillator do not form a conditional basis. 

We have also investigated the complex resonances of a typical self-adjoint operator 
by means of the standard technique of dilation analyticity. This identifies the 
resonances of the original operator with eigenvalues of any one of a family of 
associated non-self-adjoint operators indexed by an angle. The eigenvalues of these 
operators are independent of the angle, but the instability indices depend upon its 
value. 

We have discovered that for a certain operator, as is seen from Table 10, the first 
20 eigenvalues have very small imaginary parts, which is explained by the fact that 
there exists a (non-convergent) asymptotic expansion which has real coefficients 
of all orders. For higher eigenvalues the imaginary parts of the eigenvalues in- 
crease rapidly in absolute value. We have computed the instability indices of these 
eigenvalues for typical angles and discovered that they increased rapidly with the 
modulus of the eigenvalue, reaching a maximum value near the region where the 
imaginary part starts to increase (see Tables 9.1-9.3, 10). No theoretical explana- 
tion of this phenomenon exists. 

The very large size of the instability indices in both examples indicates that the 
computation of large eigenvalues of non-self-adjoint differential operators is likely 
to be intrinsically intractable in many other cases of a similar type. The same 
applies to the computation of large resonances of self-adjoint differential operators. 
The effect of rounding errors or of small perturbations of the operator may be to 
change the computed eigenvalues drastically. This discovery casts some doubt 
on the significance of theoretical investigations of the asymptotic distributions of 
resonances or of any computations of such eigenvalues for all except self-adjoint 
operators. Our experience, and that of others who work within the pseudospectral 
approach, has been that the extreme instability of large eigenvalues is the norm 
rather than a possibility which occurs only in pathological cases. 
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Figure lb. Contour map of F(z; a) 
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Figure lc. Plot and contour map of F(z; a), a = 4, b = 100 
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Figure 2a. Function F(z; 0), b = 10, 9 
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Figure 2b. Contour map of F(z; 0), b — 10, 9 — n/A 

Figure 2c. Contour map of F(z; 0), b = 10, 6 = 7r/16 
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